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ABSTRACT 

Inverse  synthetic  aperture  LADAR  (ISAL)  provides  high  resolution  images  which  are  not  limited  in  resolution  by 
the  diffraction  limit  of  the  telescope  collecting  the  imaging  data.  This  makes  it  an  ideal  method  for  imaging 
Geosynchronous  Earth  Orbit  (GEO)  space  objects.  Although  LADAR  is  a  well  researched  topic,  little  has  been 
published  on  the  suitability  of  ISAL  for  imaging  GEO  space  objects  and  the  associated  SNR  models.  This  paper 
derives  an  expression  for  the  return  signal  detected  by  a  coherent  ISAL  system.  Using  tomographic  techniques 
common  to  synthetic  aperture  radar  (SAR),  a  model  is  developed  for  the  signal-to-noise  ratio  of  the  reconstructed 
image.  Noise  sources  considered  include  photon  noise  and  laser  speckle  noise. 

1.  INTRODUCTION 

Spotlight-mode  Synthetic  Aperture  Radar  (SAR)  allows  high  resolution  terrain  imaging  with  microwaves  using 
reasonably  sized  apertures.  Similar  concepts  may  be  applied  to  laser  radar  (LADAR)  for  imaging  Geosynchronous 
Earth  Orbit  (GEO)  space  objects  using  meter  class  apertures.  Due  to  the  distance  of  GEO  objects  (around  40Mm), 
conventional  imaging  systems  require  large  apertures  to  resolve  objects.  Apertures  on  the  order  of  10m  or  more  are 
required  to  provide  useful  information.  Synthetic  aperture  systems  utilize  relative  movement  between  the  object 
and  observer  to  synthesize  larger  apertures.  For  this  reason,  much  smaller  optics  can  be  used,  thereby  reducing 
costs.  Inverse  synthetic  aperture  LADAR  (ISAL)  for  imaging  GEO  space  objects  relies  on  small,  relative 
movements  of  the  object  to  form  high  resolution  images  using  small  telescopes.  Just  as  large  telescopes  have  their 
draw  backs,  so  does  ISAL.  The  difficulty  and  costs  of  obtaining  a  10m  class  telescope  have  been  traded  for  the 
difficulty  and  cost  of  actively  illuminating  a  GEO  object  and  coherently  detecting  the  return  signal.  Therefore,  it  is 
important  to  fully  understand  the  expected  signal  and  noise  levels  for  such  a  system.  This  paper  derives  the  signal- 
to-noise  ratio  (SNR)  for  ISAL  using  first  principles  and  an  image  processing  viewpoint.  The  robust  model  can  be 
used  to  predict  future  system  performance  and  as  an  input  for  modeling  ISAL  post  processing  techniques. 


2.  BACKGROUND 

Many  imaging  techniques  using  SAR  and  LADAR  have  been  proposed.  ISAL  has  three  key  characteristics  which 
differentiate  it  from  other  techniques,  making  it  suitable  for  imaging  GEO  objects.  First,  tomographic  imaging 
techniques  are  used.  Tomography  has  typically  been  associated  with  computer-aided-tomography  (CAT)  imaging 
found  in  medical  applications  [1],  Transmission  tomography  measures  the  energy  that  passes  through  an  object  to 
obtain  a  2-D  absorption  profile  of  a  3-D  object.  Fig.  1(a)  shows  how  the  measured  absorption  creates  a  1-D 
projection  of  the  2-D  object.  This  projection  is  related  to  the  object’s  Fourier  transform  by  the  Fourier  Slice 
Theorem  [1].  It  states  that  the  Fourier  transform  of  a  1-D  projection  taken  at  angle  9  is  equal  to  a  line  (or  slice) 
through  the  Fourier  domain  of  the  2-D  object.  The  line  passes  through  the  origin  at  an  angle  9.  Therefore,  obtaining 
object  projections  for  an  infinite  number  of  angles  will  allow  perfect  reconstruction  of  the  2-D  Fourier  transform, 
and  therefore  the  2-D  object  absorption  profile  from  1-D  data.  Fortunately,  a  finite  number  of  projections  can  be 
used  for  practical  use.  Tomography  allows  3-D  images  to  be  formed  of  an  object’s  interior  non-intrusively,  making 
it  very  useful  in  medical  imaging.  Reflection  tomography  obtains  2-D  projections  of  a  3-D  object  by  measuring  the 
reflected  energy.  Fig.  1(b)  shows  how  the  energy  reflected  from  each  point  on  the  object  creates  a  projection  similar 
to  transmission  tomography.  The  differences  are  that  the  projection  obtained  by  reflection  is  from  the  cross-range 
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direction  (perpendicular  to  line-of-sight  vector),  and  the  profile  obtained  is  only  for  the  surface  of  the  object. 
However,  the  data  can  still  be  processed  the  same  way  to  reconstruct  a  3-D  object  profile  from  2-D  reflection 
profiles. 


(a)  (b) 

Fig.  1.  Transmission  tomography  (a)  measures  the  absorption  profile  of  an  object  to  recover  a  1  -D  projection. 
Reflection  tomography  measures  the  reflected  signal  as  a  function  of  range  to  obtain  a  1-D  projection. 

References  [2,  3]  show  that  tomography  can  be  applied  to  SAR.  This  allows  methods  and  algorithms  developed  for 
CAT  to  be  used  in  SAR.  In  addition,  the  tomographic  formulation  simplifies  SAR  for  the  non-radar  community  and 
allows  speckle  phenomenon  to  be  easily  explained  [22].  Tomography  has  also  been  applied  to  LADAR  [4], 

The  second  important  characteristic  that  differentiates  ISAL  from  other  LADAR  imaging  techniques  is  coherent 
illumination  and  heterodyne  detection,  instead  of  direct  detection  [5].  This  allows  the  object  to  be  reconstructed 
from  Fourier  data  which  is  offset  from  the  origin,  and  obtained  with  a  relatively  small  range  of  look  angles.  Image 
reconstruction  from  offset  Fourier  data  is  only  possible  for  complex  valued  signals  (coherent)  [2].  Consider  the 
complex  object  reflectivity  and  corresponding  Fourier  transform  given  by 

r(x,y)  =  a(x,y)exp[/0(x,y)],  (l) 

R(u,v )  =  A(u,  v)  ®  F{exp \j<p(x,y)]},  (2) 

where  a(x,y )  is  the  amplitude,  (j) (x,  y)  is  the  phase,  and  ®  indicates  a  convolution.  If  the  phase  is  highly  random 
due  to  surface  fluctuations  on  the  order  of  a  wavelength,  as  seen  in  laser  speckle,  the  Fourier  transform  of  the 
exponential  will  be  very  broad.  The  convolution  then  spreads  the  amplitude  information  across  a  wide  spectrum. 
Therefore,  offset  Fourier  data  will  contain  baseband  (origin  centered)  information  and  a  speckle  image  can  be 
reconstructed  from  this  information  [2], 

The  final  characteristic  is  that  the  ISAL  waveform  is  chirped  or  linearly  frequency-modulated  (FM).  Consider  the 
outgoing  signal 

S(t,x,y )  =  Re  { Es(x,y )  *  exp  [i(2nf0t  +  /?t2)]},  (3) 

where  Es (x,  y)  is  the  amplitude,  /0  is  the  center  frequency,  and  /?/tt  is  the  chirp  rate  in  Hz/s.  The  frequency,  or  first 
derivative  of  the  phase,  changes  linearly  with  time.  There  are  two  advantages  to  using  a  chirped  waveform.  First,  it 
has  a  large  time-bandwidth  product.  The  bandwidth  is  given  by 

Bc  =  Ptc/tz,  (4) 

where  rc  is  the  chirp  duration  [3],  From  Eq.  (4),  it  is  clear  that  long  pulses  result  in  high  bandwidths.  This  is  an 
opposite  result  when  compared  to  a  continuous  wave  (CW)  short  burst,  where  a  short  pulse  duration  leads  to  high 
bandwidth.  The  large  time-bandwidth  product  means  that  more  energy  can  be  transmitted,  resulting  in  higher  SNRs. 
The  other  advantage  of  a  chirped  waveform  is  that  the  processed  return  signal  will  be  the  Fourier  transform  of  the 
object’s  reflectivity  projection,  offset  from  the  origin  as  shown  in  Fig.  2.  Image  reconstruction  will  be  possible  from 
this  offset  data  due  to  the  random  phase.  The  benefit  of  collecting  offset  Fourier  data  is  apparent  from  the  figure.  A 


more  complete  UV  fill  is  obtained  from  the  offset  data,  when  compared  to  the  baseband  data.  Baseband  SAR  and 
SAL  systems  require  a  wide  range  of  look  angles  to  reconstruct  a  useful  image.  For  Low  Earth  Orbit  (LEO)  space 
objects.  Ad  will  be  large,  and  baseband  systems  can  be  used  [4],  However,  the  relative  changes  in  angle  for  GEO 
objects  are  small.  For  this  reason,  ISAL  is  well  suited  for  GEO  imaging. 


Fig.  2.  Spatial  frequency  UV  fill  for  look-angles  in  the  range  Ad.  Only  a  small  cone  around  the  origin  gets  filled  in 
for  the  baseband  data.  The  offset  data  obtained  by  ISAL  is  more  complete. 


3.  ANALYSIS 

The  following  analysis  is  based  on  spotlight-mode  SAR  tomographic  concepts  presented  in  [3].  A  similar  analysis 
has  been  applied  to  SAL  for  planet  imaging  from  satellite  systems  [8],  Although  Lucke  and  Rickard  help  bridge  the 
gap  between  the  radar  and  optical  imaging  communities,  their  analysis  is  simplified  and  therefore  their  models  lack 
generality  and  robustness  [8],  This  analysis  seeks  to  develop  a  SNR  model  of  ISAL  in  a  general  manner,  using  first 
principles  where  possible. 


3.1  Basic  Formulation 


Consider  the  2-D  transmitter/receiver  geometry  shown  in  Fig.  3  with  a  transmitted  signal  given  by  Eq.  (3).  The 
instantaneous  signal  at  the  receiving  aperture  reflected  from  points  between  —z1  and  z1  on  the  object,  is  given  by 

Srefiectedit.y )  =  Re  [ a  *  o (y, z)  *  Es[t  -  t(z), y]  *  exp  [t2?r  [/0(t  -  r(z))  +  /?(t  -  r(z))2]}  dz],  (5) 


where  a  is  a  scaling  coefficient  accounting  for  signal  loss  due  to  atmospheric  absorption,  beam  divergence  to  and 
from  the  object,  and  any  transmission  loss  due  to  optics.  o(y,  z)  is  the  object  field-reflectivity  density  (reflectivity 
over  a  small  incremental  range  Sz,  centered  at  z),  r(z)  =  is  the  propagation  time  from  point  z  to  z  =  0.  The  round 
trip  propagation  time  to  the  center  of  the  object  and  back,  given  by  r0  =  2  z0/c,  has  been  left  off  for  convenience.  It 
is  important  to  note  that  the  chirp  duration  rc  is  much  longer  than  the  round  trip  propagation  delay  across  the  object 
and  back  given  by  rp  =  2 zx/c,  where  c  is  the  speed  of  light.  Equation  (5)  is  only  valid  during  the  time  [3] 


+  <t<  tJL-X 

2  2  2  2 


(6) 


During  this  period,  the  return  signal  contains  information  from  the  entire  depth  of  the  object  and  can  be  considered  a 
convolution  between  the  object’s  reflectivity  profile  and  the  transmitted  pulse. 


y 


Fig.  3.  Geometry  used  in  ISAL  analysis,  along  with  chirped  waveform  over  duration  rc. 

For  heterodyne  detection,  the  return  signal  is  split  and  mixed  with  a  local  oscillator  (LO)  which  consists  of  in-phase 
and  quadrature  versions  of  the  FM  chirp,  given  by 

C,(t)  =  ^Elo  cos[27r/0(t  -  Td)  +  Pit  -  rd)2],  (7) 

CQ(t)  =  -  \el0  s\n[2nf0(t  -  Td)  +  P(t  -  Td)2\  (8) 

where  Td  is  a  deliberate,  relative  time  delay  between  the  LO  and  returning  signal  and  Elo  is  the  amplitude  of  the  LO. 
The  time  delay  allows  the  temporal  frequency  content  of  the  mixed  fderamped)  signal  to  be  controlled.  The  Vi  in 
Eqs.  (7)  and  (8)  result  from  one  LO  being  split  between  the  two  detectors  (in-phase  and  quadrature). 


3.2  Analysis  of  In-Phase  Intensity 


The  following  section  considers  the  in-phase  channel.  A  similar  analysis  can  be  applied  to  the  quadrature  channel. 
After  mixing,  the  intensity  on  the  detector  is  given  by 


+2 


hit, y)  =  | \sreflected{t,y)  +  C,(t)|  , 

h(f,y)  =  \^elo  *  cos[2nf0{t  -  rd)  +  Pit  -  rd)2]j 

1  fZl 

-  a  *  o(y, z)  *  Es  (t,y)  *  cos(27r{/0[t  -  t(z)]  +  p[t  -  r(z)]2})  dz 
Z  ■'“^1 

rzi  \ 

J  a  *  oiy,z )  *  -JmiEL0  *  Es  (t,y)  *  cos[27r/0(t  -  td)  +  pit  -  rd)2]  *  cos(27r{/0[t  -  r(z)]  +  p[t  -  r(z)]2})  dz 


(9) 


(10) 


where  rjh  is  the  heterodyne  mixing  efficiency  [8].  The  first  term  in  Eq.  (10)  results  in  a  direct  current  (DC)  term  and 
a  cosine  with  frequencies  centered  at  2<u0,  which  our  detector  will  not  respond  to.  Therefore,  the  first  term  will  be 


hit.y) 


1st  term 


(ii) 


The  second  can  be  written  as 


hit, y)2nd  term  =  [g*gsO->9]  iy ,  z)2  dz  +  o(y,  z)o  (y,  z)  *  cos[(z  -  z)^i(0  o  -  2^t)]  dzdz],  (12) 

where  z  is  a  dummy  variable  used  for  integration.  At  any  time  during  the  processing  window,  the  second  term  gives 
a  DC  component  for  each  point  on  the  target  weighted  by  the  reflectivity  squared.  It  also  contains  cross  terms 
between  all  points  on  the  object,  with  a  frequency  which  depends  on  the  separation  of  the  two  respective  points. 

The  beat  frequencies  which  occur  in  this  term  range  from  zero  where  z  =  z,  to  ±  Qfiz1/c,  when  the  reflected  signal 


from  the  far  ends  of  the  object,  ±z1  interfere.  As  will  be  seen  later,  this  is  the  same  frequency  range  as  the  third 
term  which  contains  the  desired  signal.  For  GEO  imaging,  the  return  signal  will  be  very  weak  compared  to  the  LO 
(a  *  o(y,z)  *  Es{t,y )  «  Ei0)  and  the  first  and  third  terms  of  Eq.  (10)  will  dominate.  For  this  reason,  the  second  term 
can  be  ignored.  In  addition,  the  time  delay  xd  can  be  used  to  shift  the  frequencies  of  the  desired  signal  so  that  the 
two  terms  don’t  overlap  in  bandwidth. 

Finally,  we  consider  the  last  term  of  Eq.  (10).  The  two  cosines  can  be  expanded  into  two  terms  with  frequencies 
given  by 

fA=fi-f2=^-2Prd,  (13) 

fB=fi+f2  =  2(o0+4pt-2pTd-^.  (14) 

The  lowest  frequency  will  occur  at  the  beginning  of  the  processing  window,  given  in  Eq.  (6).  Substituting  the  lower 
value  into  Eq.  (14)  gives 

fB  =  2o)0-2/?[Tc-(rp-Td)]-^.  (15) 

Assume  that  the  chirp  duration  will  be  much  longer  than  the  round  trip  propagation  time  across  the  object  and  the 
deliberate  time  delay  combined.  That  is,  xc  »  (xp  —  rd).  The  second  term  in  Eq.  (14)  becomes  2nBc.  Since  the 
bandwidth  is  much  smaller  than  the  center  frequency,  this  term  is  negligible  compared  to  the  first.  The  last  term  will 
also  be  negligible  compared  to  the  first.  This  leaves 


/b  ~  2 cog.  (16) 

Once  again,  our  detector  will  not  respond  to  such  high  frequencies.  Therefore,  only  the  cosine  term  with  frequencies 
in  the  range  of  fA  will  be  detected.  The  third  term  of  Eq.  (10)  reduces  to 

h(t,y)3rd  term  =  «  *  o(y, z)Jr^EL0  *  Es  (t,y)  *  cos ((z  -  ^Td)  ^  [oj0  +  2 /3t]  -  +  /?rd2}  dz,  (17) 

The  quadratic  phase  term  4/?z2/ c2  acts  to  blur  the  image  in  the  cross  range  dimension  [3],  Through  an  example, 
Jakowatz  states  that  as  long  as  this  term  is  below  7r/4,  the  blurring  can  be  neglected.  Even  for  relatively  large  GEO 
objects  (50m)  and  high  chirp  rates  (3GHz/s),  this  term  will  only  be  a  few  mrad  and  can  be  ignored.  Finally,  the  last 
term  in  the  phase  can  also  be  neglected  if  it  is  7r/4  or  less.  This  means  the  time  delay  is  limited  to 


1^1  <J%- 


(18) 


Eq.  (17)  then  becomes 

h(t,y) srdterm  =  a^EL°Es{t'y)  0(y,  z)  *  COS  [(z  ~  ^d)  \  No  +  2/?t]}dz, 


(19) 


This  term  can  be  thought  of  as  the  real  part  of  the  Fourier  transform  of  the  object  reflectivity  with  an  additional 
phase  term  which  is  determined  by  xd.  If  we  pull  this  additional  phase  term  out  of  the  integral  over  z,  we  can  write 


y)3 


aJnhEL0Es(.t,y) 


Re  | exp  [/zZ]  o(y,z)  *  exp  (— /zZ)  dzj. 


(20) 


I,(t,yhrdterm  =  a^hELfsit'y) Re{exV[- jzZ]  0(y,Z)}, 


(21) 


C  2 

where  z  =  -rd,  and  Z  —  -  (e>0  +  2 fit).  0(y,  Z)  is  the  Fourier  transform  of  the  object’s  reflectivity-density  function, 
for  a  given  y  [3],  The  spatial  frequency  Z  is  a  function  of  time.  Therefore,  at  each  time  t  during  the  processing 
window,  the  signal  corresponds  to  the  object’s  Fourier  transform  evaluated  at  a  particular  spatial  frequency.  Over 
the  processing  window,  the  range  of  Z  is  given  by 


(22) 


l  [«j0  -  PTc\  <  z  <  -c  [tO0  +  Z?TC], 


and  therefore  the  spatial  frequency  bandwidth  is 


AZ  = 


(23) 


This  offset  Fourier  data  is  centered  at  2 (j)0/c. 

3.3  Detector  Output 

The  current  leaving  the  in-phase  detector  is  given  by  [8] 

i/(0  = 

where  r/d  is  the  detector  quantum  efficiency,  qe  is  the  charge  of  an  electron,  h  is  Planck’s  constant,  and  v  is  the 
optical  frequency.  If  the  second  term  is  ignored  in  Eq.  (10),  the  current  becomes 

i/(0  =  f^fJ"  [lf  +  \S-Zla  *  o^^Elo  *  Es  ( t,y )  *  cos[(z  -  jTd)^[w0  +  2/?t]j  dz]  dy. 

For  simplicity  assume  that  the  amplitude  of  the  return  signal  is  uniform  across  y,  and  thus  the  integrated  current 
becomes 

ii(t)  =  | j^Adet  l^f-  +  lf-Zl  a  *  o(z)4dhELO  *  Es  (t)  *  cos  [(z  -C-id)-c  [w0  +  2/?t]}  dz], 

where  Adet  is  the  area  of  the  detector.  The  average  electric  field  amplitudes  are  related  to  the  average  number  of 
photons  N  by 


(24) 


(25) 


(26) 


E  = 


2hvN 

Adetz 


(27) 


where  the  integration  time  is  r  =  ( rc  —  rp)  [8].  Using  Eq.  (27),  the  current  becomes 

h(t)  =U^[nlo  +  2a*Jrj ^Vlo^s  *  *  cos  j(z  -  z')  +  2j?t]}dz],  (28) 

Finally,  the  current  from  the  in-phase  and  quadrature  detectors  can  be  combined  to  represent  a  complex  signal  given 
by 

i(t)  =  ^[NloQ-  +;)  +  2a  *  exp[— ;(z  +  z')^[co0 +  2)?t]}dz],  (29) 


3.4  Processing  the  Detector  Output 

After  the  bias  is  subtracted,  the  signal  represents  a  scaled  and  phase  shifted  Fourier  transform  of  the  object’s 
reflectivity  density  function.  To  recover  the  object’s  reflectivity  profile,  and  therefore  the  image,  the  inverse  Fourier 
transform  is  taken.  Equation  (29)  can  be  written  to  group  the  time  varying  factors  in  the  phase 


i(t)  =  k  o(z)  *  exp  (—j2n  t  +  dz, 

where  l  is  the  mean  removed  current,  /?  =  [1  /2n,  and 

k  =  a-^T^n*NL0*Ns. 


(30) 


(31) 


A  variable  substitution  can  be  made  where 


f  = 


4/?(z+z) 

c 


(32) 


is  the  temporal  beat  frequency  from  point  z  on  the  object,  and  Eq.  (30)  becomes 

A<0  =  k  f-/i  *  exP  [~P 2n  (fl  +  )]  jfdf- 

o(z)  is  the  reflectivity  per  Sz  increment,  where  Sz  —  zx  —  z2.  When  changing  variables,  this  becomes  the 
reflectivity  per  Sf  given  by 

<A  =  0(A-/2). 


(33) 


(34) 


The  c/4/?  factors  in  the  integral  cancel  each  other  out.  We  now  take  the  inverse  Fourier  transform  the  signal  given 

by 

1  (/)  =  k  4  exP  [~i2n  +  70\  d/exP (2 njft)dt. 


where  /  is  a  dummy  variable  used  for  integration.  Now  the  order  of  integration  is  switched  and  the  exponentials  are 
combined  to  get 


A/)  =  k  4  o  (g  -  z)  14  exp  [-j2n  (ft  +  g)]  dtdf. 


(36) 


Next,  the  non-time-varying  terms  in  the  exponential  is  moved  outside  the  integral. 

A/)  =  k  4  o  (0  -  z)  exp  (a'Ztt AAl)  J_4  exp(-j2n(f  +  f)t )  dtdf , 

The  integration  can  now  be  carried  out. 

A  f)=k  /A  o  (0  -  z)  exp  (i27T0)  42{cos[2tt(/  +  /)t]  +  t  sin[27r(/  +  /)t]}dtd/. 

For  convenience,  we  can  multiply  the  last  factor  by  t/t  to  create  a  sine  function 

l (/)  =  Tfc  4  0  (0  -  x)  exP  (i2n 0)  sinc[7rr(/  +  f)\df. 


(37) 


(38) 

(39) 


(40) 


Had  the  Fourier  transform  been  taken  for  an  infinite  signal  (  i.e.  r  =  oo),  the  sine  function  would  become  a  delta 
function  with  zero  width.  This  delta  function  would  perfectly  sample  the  object  reflectivity  for  each  frequency  over 
the  given  integral.  However,  the  sine  has  a  central  lobe  width  of  2/t  with  "ringing”  past  that.  This  acts  to  blur  and 
distort  the  sampled  frequencies.  Since  these  frequencies  represent  different  ranges  on  the  object,  the  reconstructed 
image  is  blurred  in  range.  Using  the  derivative  of  Eq.  (32)  we  can  relate  the  blur  in  frequency  to  the  blur  in  range  by 
substituting  dfsinc  —  2/r. 


dz 


sine 


(2 /t)c 
4(/?/2tt) 


cn 

r/?’ 


(41) 


Figure  4  shows  how  the  sine  function’s  central-lobe  width  changes  with  integration  time  for  two  different  chirp 
rates.  The  blur  in  range  decreases  as  the  integration  time  increases.  However,  it  is  important  to  note  that  the 
resolution  will  also  increase.  Reference  [3]  gives  the  range  resolution  as 


(42) 


If  we  divide  the  sine’s  width  by  the  range  resolution,  the  result  is 


dzsinc  _  2tc  (43) 

Sz  t  ’ 


where  Eq.  (4)  has  been  used.  Finally,  rc/ (rc  —  rp)  »  rc,  thus  the  blur  caused  by  the  finite  integration  time  will 
always  result  in  a  blur  of  two  pixels. 


Integration  time  [s] 


Fig.  4.  Blur  in  range  z  as  a  function  of  integration  time  for  two  different  chirp  rates. 


As  the  integration  time  increases,  the  sine’s  width  will  be  small  compared  to  significant  changes  on  the  object’s 
profile.  Therefore,  the  sine  function  can  be  modeled  as  a  Dirac  delta  function  and  Eq.  (40)  becomes 


I  (/)  =  kT  J_A  o  (0  +  z)  exp  [jin  ^jf)8(f  +  f)df. 


(44) 


Carrying  out  the  integration  leaves 


/(/)  =  kT  *  o  (g  +  z)  exp 


(45) 


where  the  integration  over  /  has  changed  o(/)  from  a  reflectivity  density,  to  a  unit-less  reflectivity.  Equation  (32) 
is  used  to  change  variables  back  to  z,  giving 


/(z)  = 


a*ridqe 


y/rih*NLO*Nso(z)  exp  (;  [2(z+^°]}, 


(46) 


The  difference  between  /(/)  and  /(z)  can  be  confusing.  The  two  are  the  same  function,  just  evaluated  at  different 
frequencies.  7(z)  is  simply  1  (j  =  —z'j.  The  end  result  is  a  scaled  and  phase-shifted  version  of  the  object 
reflectivity  along  the  z  axis.  The  last  factor  applies  a  phase  to  each  range  element.  This  results  from  the  time  delay 
and  shifting  the  signal  to  baseband  prior  to  the  inverse  Fourier  transform.  Finally,  we  estimate  1-D  object 
reflectivity.  To  do  this,  the  extra  factors  are  divided  out  and  the  magnitude  squared  is  taken. 


o(z)  = 


_ Tf) _ 

(^)*{VdqeJrih*NLo) 


(47) 


o(z)  =  Nso2(z). 


(48) 


Here  we  see  that  the  signal  of  interest  is  simply  the  energy  spectral  density  (ESD)  of  the  heterodyne  output.  This 
will  be  important  when  developing  the  stochastic  model. 


3.5  Reconstructed  Image  Using  Backprojections 


Although  o(z)  is  a  one -dimensional  function  of  range,  it  actually  represents  a  projection  of  the  two  dimensional 
object  reflectivity  function  taken  at  angle  0  [3],  A  subscript  dg  (z)  is  added  to  note  the  angle  for  which  the  return 
signal  was  measured.  Figures  5  and  6  illustrate  this  concept.  For  angle  6 ,  the  LOS  vector  from  the  aperture  to  the 
object  becomes  z’.  This  is  the  direction  of  illumination  and  the  slice  of  Fourier  data  contained  in  the  return  signal 
shown  in  Fig  5(c)  is  only  a  function  ofz’.  Figure  6  provides  the  1-D  perspective  of  the  object,  return  signal,  and 
processed  signal  shown  in  Fig.  5.  Figure  5  shows  how  this  1-D  return  signal  is  actually  a  line  in  the  2-D  spatial 
frequency  domain.  In  Fig.  5,  the  2-D  inverse  Fourier  transform,  or  backprojection  [1]  representing  the  processed 
signal  for  look  angle  6,  shows  how  the  1  -D  processed  signal  is  “smeared”  across  the  y-z  plane  at  an  angle  6.  All 
points  along  lines  parallel  to  y’  have  the  same  value  for  a  single  look  angle. 


Fig.  5.  (a)  2-D  circular  object  in  y-z  plane  with  rotated  axes  z’  and  y’  shown,  (b)  Fourier  transform  of  (a),  (c)  single 
slice  of  (b)  corresponding  to  the  return  signal  along  the  z’  axis,  and  (d)  inverse  Fourier  transform  of  (c)  with  rotated 

axes  shown. 
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Fig.  6.  (left)  Simple  1-D  example  of  rectangular  object  (slice  of  2-D  circle  in  Fig. 5),  (middle)  Fourier  transform  of 

object,  and  (right)  inverse  Fourier  transform  of  (middle) 

For  each  look  angle,  a  different  backprojection  is  obtained.  When  combined,  more  information  about  the  image  is 
gained.  Figure  7  demonstrates  how  the  summation  of  backprojections  (bottom  row)  improves  the  image.  From 
these  simple  examples,  it  is  clear  that  each  pixel  value  in  the  z-y  image  plane  is  the  sum  of  individual 
backprojections.  We  must  mathematically  determine  the  total  contribution  to  each  pixel  in  the  y-z  plane.  First 
consider  a  single  return  and  backprojection.  For  a  given  1  -D  return  we  express  the  value  of  pixel  (y,z)  in  the  final 
image.  The  1-D  projection  is  smeared  along  the  6+n/2  direction  (orthogonal  to  z’  axis)  [1],  The  value  of  the  1-D 
return  at  point  z’  is  given  to  each  pixel  in  the  2-D  array  which  lies  along  the  line  z  =  (tan  6)y  +  (sec  0)z'.  For  a 
pixel  at  location  (y1(  z±)  in  the  2-D  image  \g(y,  z)  |2,  the  value  contributed  from  a  single  backprojection  is  given  by 

\ge(.yi'Zi)\2  =  I og(z'  =  z1cosB  -yiSintf)!2  =  |A/riV^o(z1  cos 6  -y1  sin 0 ) exp  [j2(Zl c°s 9~*1  sm 6)a)°] |  .  (49) 

The  value  of  each  pixel  in  the  final  reconstructed  image  (composed  of  all  measured  look  angles)  is  given  by 


2 


(50) 


|  o (y,  z) | 2  =  IE"®  -  JTso(zi 


fi2fzi  cos  6- 


where  the  angle  0  —  rig  *  SB,  Ng  is  the  total  number  of  look  angles,  and  56  —  Ng/rig.  For  additional  information 
about  techniques  such  as  convolution  backprojection  (CBP)  which  can  improve  the  final  image  quality,  the  reader 
may  consult  [1,  3,  9], 


Fig.  7.  From  the  object  in  Fig.  5,  slices  of  the  2-D  spatial-frequency  domain  for  different  look  angles  (top  row)  are 
inverse  Fourier  transformed  (middle  row),  and  added  together  to  form  an  image  (bottom  row). 

3.6  Signal-to-Noise  Ratio 

The  dominate  noise  sources  for  heterodyne  detection  are  laser  speckle  noise  and  shot  noise.  Both  background  and 
detector  noise  can  be  neglected  [10].  As  will  be  shown,  the  electro-magnetic  (EM)  field  that  is  incident  on  the 
detector  has  circularly  complex  Gaussian  statistics  and  the  detector  output  has  Poisson  statistics  which  depend  on 
the  speckle  field  intensity  [11,  12].  The  detector  output  can  therefore  be  modeled  as  a  doubly  stochastic  Poisson 
process  as  described  by  [12].  A  doubly  stochastic  Poisson  process  is  a  Poisson  process  with  a  random  rate  function. 
For  ISAL,  the  Poisson  rate  function  A(t)  is  the  complex  detector  output  given  by  Eq.  (29),  divided  by  qe  to  get  the 
number  of  photo-electron  events. 

_  h(t)+jiQ(t)  (51) 

Qe 

If  the  z  and  z'  terms  in  the  cosine  and  sine  of  Eq.  (51)  are  expanded,  the  z’  terms  can  be  moved  out  of  the  integral, 
and  the  rate  function  can  be  written  as 

A(t)  =  +;)  +  *  NL0  *  /V,exp  (w0  +  2/?t)]  o{z)  {cos  (w0  +  2 /St)]  +  j  sin  (<u0  +  2 /St)]}  dz],  (52) 

If  we  consider  integrating  over  a  range  dz  large  compared  to  a  wavelength  but  small  compared  to  changes  in  the 
signal  envelope,  we  can  model  the  return  signal  as  a  planar-surface  laser  speckle  problem  [13]  which  is  thoroughly 
treated  by  [11,  12],  The  random  fluctuations  in  z  of  the  rough  object  surface  will  create  randomness  in  the  return 

signal.  In  Eq.  (52)  both  o(z)c os  (a>0  +  and  o(z)  sin  (a>0  +  represent  random phasor  sums 

integrated  over  y  in  Eq.  (26),  and  can  be  considered  zero-mean  independent  and  identically  distributed  Gaussian 
random  variables  [12],  assuming  that  2/?t  «  a>0.  The  complex  rate  function  is  therefore  a  circular  complex 
Gaussian  random  variable,  with  mean  JVL0(  1  +  j). 


The  rate  function  and  corresponding  ESD,  EA(/)  can  be  expressed  as 


^(t)  —  ~Nlo(  1  +  7)  +  i(t)  , 

oT 

Ea(/)  =  E{|A(/)|2}  =  E{|^iVL0(l  +  ;)«5(/)  + 

Ea(/)  =  E{A(/)A*(0)  =  E{[^AiL0(l  +;)5(/)  +  /(/)]  [^NL0(  1  - 7)5(0  +  HO]}, 

Ea(/)  =  e{2  (^Wlo)2  5(/)5(0  +  ^JVl0(1  -7)5(0/ (/)  +^^lo(1  +;)5(/)r(0  +/(/)/* (o}, 
Ea (/)  =  2  (^Nlo)2  S(f)S(0  +  ^~tNlo(1  -;){5(0E[/(/)]  +  5(/)E{/‘(0)}  +  E{/(/)/‘(0), 


where 


E{7 (/)}  =  E  *  Nlo  *  Nso  (0  -  z)  exp  (i2?r  g)], 


(53) 

(54) 

(55) 

(56) 

(57) 

(58) 


Just  as  we  did  above  to  model  the  laser  speckle  in  Eq.  (52),  we  can  expand  Eq.  (58)  in  a  similar  way  using  Eq.  (32). 
Therefore,  Eq.  (58)  becomes 


W (/))  =  ^  yjrin  *  Nlo  *  Ns  exp  (gg)  E  (o  (z)  exp  (gg)}, 


(59) 


Just  as  before,  the  factors  inside  the  expectation  can  be  treated  as  a  circular  complex  Gaussian  random  variable,  and 
E {/(/)}  =  0.  Next,  let  us  find  the  second  moment. 


where 


E{/(/)r(0)  =  E  {t2o  (fif  -  z)  exp  exp  (-1  ^)J 

k=^h.  p(!Hf). 

For  /  4=  the  second  moment  is 


E{/(/)/‘(0)  =  fe2E {o  g  +  i) exp  (-i §)} E(og  +  z) exp (i ^)}  =  0. 


When  /  =  (, 


E{/(/)/*(0)  =  k2o2  g  +  z)  S(/  -  0- 


Finally  Eq.  (57)  becomes 


E 


i(/)  =  2  (0A1lo)2  5(/)  +  k2o2{z). 


(60) 

(61) 

(62) 

(63) 

(64) 


Now  that  the  rate  function  has  been  analyzed,  we  can  consider  shot  noise.  The  ESD  of  the  doubly  stochastic 
Poisson  process  will  give  the  mean  signal  in  the  image  for  a  single  projection  (prior  to  the  bias  terms  being 
removed).  In  this  case  the  output  of  each  detector,  in-phase  and  quadrature,  represents  a  complex  Poisson  impulse 
response  with  the  complex  signal  expressed  as  [12] 


and  Fourier  transform 


i(0  =  En=i5(t-tn)+7  2fc=i5(t-tfc). 


7(/)  =  Hn=i  exp(j2nf  tn)  +  )l£=iexp(/27r/tfc). 


(65) 

(66) 


Goodman  [12]  shows  that  the  energy  spectral  density  for  a  doubly  Poisson  process  can  be  expressed  as 

E  3(0  =  %+  Ea  ,(/)  +  Kq  +  EAq(/),  (67) 


where  KI  and  Kq  are  the  average  number  of  photo-electrons  generated  by  the  in-phase  and  quadrature  stochastic 
rate  functions,  respectively.  If  we  assume  identical  detectors,  then  Eq.  (67)  can  be  written  as 


E,(/)  =  2K  +  EA(/).  (68) 

K  can  be  found  by  time  averaging  either  the  real  or  imaginary  part  of  the  mean  rate  function  [12] 

K  =  j^2  E[Re{/l(t)}]dt  =  E[ImU(t)}]dt  =  Re{X}r  =  Im{I}r.  (69) 

K  =  '-fNL0.  (70) 


Finally,  the  ESD  can  be  expressed  as 


EjO)  =  1JfNL0  +  2  (^Nlo)2  Sif)  +  k2o2  +  z),  (7 

This  result  shows  that  the  mean  pixel  value  for  each  frequency  (and  respective  range),  is  equal  to  the  scaled  return 
signal  weighted  by  the  object  reflectivity,  plus  a  bias  term,  representing  the  LO,  and  a  spike  at  the  origin  due  to  the 
bias  of  the  rate  function.  In  the  end,  we  wish  to  estimate  just  the  signal  of  interest.  This  is  given  by 


O(z)  = 


Vd„ 

~r nlo " 


Fr4V?7fc*iVio] 


=  (Vso2(z), 


(72) 


Now  that  the  first  moment  Ej(/)  has  been  found,  the  second  moment  of  |J(/)|2  is  needed  to  find  the  variance,  and 
corresponding  SNR.  In  this  case,  that  means  finding  the  fourth  order  moment  of  the  Fourier  transform  of  a  doubly 
stochastic  Poisson  process  expressed  as 


E{|7(/)|4}  =  E{Zn=i'Zm=i'Lp=i'Zq=iexp(j2n:ftn)exp(j2nftm)exp(j2nfp)exp(j2nftq)}.  (73) 


Fortunately  for  us.  Ref.  [12]  has  an  entire  appendix  (Appendix  C)  devoted  to  this  lengthy  calculation.  The  final 
result  is  given  by 

E{|7(/)|4}  =  K  +  2K2  +  4(1  +  /T)Ea(/)  +  Ea(2/)  +  2EA2(/),  (74) 

which  makes  the  variance 

^  2=K  +  K2  +  2KEA(f)  +  4Ea(/)  +  Ea(2/)  +  EA2(/).  (75) 


From  above,  K  for  a  single  Poisson  process  will  be  2 K  for  a  complex  Poisson  process  when  the  real  and  imaginary 
parts  have  equal  means.  Equation  (75)  becomes 

^(n\2=2  K  +  4  K2  +  4/fEA(/)  +  4Ea(/)  +  Ea(2/)  +  EA2(/).  (76) 


Here  we  see  what  Ref.  [12]  calls  the  "half-frequency  phenomenon”  in  the  fifth  term.  The  noise  at  frequency/has 
contributions  from  the  signal  present  at  2 f.  However,  this  term  can  be  ignored  if  the  beat  frequency  bandwidth  is 
smaller  than  the  center  beat  frequency.  In  this  case,  the  bandwidth  will  not  overlap  with  any  noise  at  the  half¬ 
frequencies.  For  this  to  be  the  case,  the  negative  time  delay-distance  z  must  be  much  greater  than  the  object  width. 
That  is. 


A/  <  \f0, 

z  >  3z1. 


(77) 

(78) 


Once  again,  this  shows  that  the  LO  time  delay  will  reduce  noise.  However,  the  value  of  the  time  delay  is  bounded 
by  both  Eq.  (18)  and  (78).  In  addition,  the  noise  terms  with  8 (/)  can  be  ignored  since  the  beat  frequency  /  will  not 
be  zero.  Moving  on,  we  substitute  K  and  EA(/)  into  Eq.  (76)  giving 


a\Kz)\2  =2(jf-NL°)  +  4:(ciNL°)  +4(YWLo)fc2o2(z) +  4fe2o2(z) +  fe4o4(z). 


(79) 


Since  we  divided  the  signal  estimate  by  a  constant  in  Eq.  (72),  the  variance  must  be  divided  by  the  square  of  the 
same  constant,  giving 


'P(z)l2 


_  2(^-Wio)+4(^Wio)2+4(^Wio)fc202(z)+4fc202(z)+fc404(z) 


(80) 


Finally,  the  SNR  for  a  single  projection  can  be  expressed  as 


SNRg  = 


SNRg  = 


Nso2(z ) 


2  (^-Wi0)+4(^JVi0)  +4(^Nio)fc202(z)+4fc202(z)  +  k404(z) 


(81) 

(82) 


When  the  individual  projections  are  “smeared”  across  the  final  image  in  the  backprojection  process,  the  statistics 
will  not  change.  The  above  SNR  is  applicable  to  each  pixel  in  the  final  image. 

4.  RESULTS 

It  is  now  useful  to  analyze  Eq.  (82)  as  a  function  of  signal  photons.  The  solid  blue  line  in  Fig.  (8)  was  computed 
from  Eq.  (82)  for  an  ideal  case  when  there  are  no  losses  in  the  system  ( rid  —  r/h  =  a  —  o  —  1).  The  curve  shows 
that  for  a  low  number  of  signal  photons,  the  shot  noise  dominates.  When  the  signal  increases,  the  SNR  levels  off  at 
1.  This  is  consistent  with  [11]  and  is  a  result  of  the  laser  speckle.  No  matter  how  many  signal  photons  are  received, 
the  image  SNR  will  not  increase  above  1.  There  are  techniques  to  improve  this  limit,  but  those  will  not  be  discussed 
at  a  this  time.  For  more  information,  see  [8,  3],  The  dashed  blue  line  in  Fig.  (8)  shows  how  under  less  ideal 
circumstances,  ( Tjd  =  0.4,  rih  =  0.5,  a  —  0.3,  o  =  V0.4),  the  SNR  shifts  to  the  right,  resulting  in  more  signal 
photons  being  required  to  achieve  the  same  SNR. 

Next,  we  wish  to  compare  the  above  results  with  those  obtained  in  [8],  The  green  lines  in  Fig.  (8)  show  Lucke’s 
(RL)  SNR  for  the  ideal  and  non-ideal  cases.  Both  SNR  models  have  similar  shapes,  although  the  RL  SNR  is  a  bit 
more  optimistic  than  the  one  developed  in  this  paper.  There  are  two  reasons  for  this.  First,  the  RL  SNR  does  not 
account  for  the  fact  that  the  beams  must  be  split  for  the  two  detectors.  Second,  the  RL  model  does  not  include 
object  reflectivity  or  transmission  losses  along  the  optical  path.  When  those  two  issues  are  ignored,  the  SNR  model 
given  in  Eq.  (82)  matches  the  RL  model.  The  RL  model  could  be  easily  modified  to  accommodate  these  parameters. 
The  differences  are  not  major,  but  may  give  unrealistic  results  when  not  considered. 


Comparison  of  Signal-to-Noise  Ratio  Estimates 


Fig.  8.  SNR  computed  from  Eq.  (82),  as  a  function  of  signal  photons  (CP,  blue).  SNR  computed  from  [8]  (RL, 

green).  Both  ideal  and  non-ideal  conditions  are  shown. 


Finally,  it’s  helpful  to  apply  the  SNR  model  described  in  Eq.  (82)  to  see  how  ISAL  sample  images  might  appear  for 
different  signal  levels.  Figure  9  shows  eight  images  over  a  range  of  SNR  levels.  A  no-noise  image  was  formed  by 
taking  100  radial  slices  through  the  Fourier  transform  of  the  pristine  image.  Then,  a  convolution  backprojection 
algorithm  [1]  was  used  to  reconstruct  a  no-noise  ISAL  image.  Finally,  Gaussian  noise  was  added  to  each  pixel  with 
a  mean  given  by  the  corresponding  signal  level  Nso2  (z),  and  the  variance  was  computed  from  Eq.  (82).  A  Gaussian 
distribution  was  selected  based  on  [8].  The  figure  highlights  how  ISAL  images  differ  from  conventional  images. 

As  the  SNR  increases,  the  resolution  stays  constant;  however  the  ability  to  distinguish  the  object  improves. 
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Fig.  9.  Example  image  of  C  SAT  for  various  SNR  levels. 


5.  CONCLUSIONS 


Through  the  above  analysis,  a  SNR  model  for  ISAL  images  has  been  developed  using  terminology  and  techniques 
common  to  the  image  processing  community.  Tomographic  methods  were  used  to  model  and  explain  how  chirped 
coherent  illumination  results  in  the  detection  of  offset  Fourier  data  in  the  spatial  frequency  domain.  Since  the  laser 
speckle  acts  to  modulate  baseband  data  across  the  Fourier  domain,  it  allows  us  to  reconstruct  the  object  profile  from 
offset  data.  This  aspect  of  ISAL  differentiates  it  from  other  LADAR  imaging  techniques.  A  small  range  of  look 
angles  can  be  used  to  image  GEO  space  objects.  While  laser  speckle  helps  us  in  that  sense,  it  also  limits  the 
achievable  SNR  to  one.  In  low  signal  levels,  we  saw  that  shot  noise  dominates.  As  the  signal  increases,  the  speckle 
dominates,  and  the  SNR  levels  off  at  unity.  Methods  to  improve  ISAL  SNR,  will  be  the  subject  of  future  research. 
The  model  developed  in  this  paper  will  provide  a  useful  tool  when  designing  an  ISAL  system  for  GEO  imaging  or 
when  modeling  other  ISAL  aspects  such  as  transmitter  power,  or  image  improvement  algorithms. 
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